###########################
### IMPERIAL RULE, THE IMPOSITION OF BUREAUCRATIC INSTITUTIONS, AND THEIR LONG-TERM LEGACIES
### WORLD POLITICS, VOL. 71, NO. 4 (OCTOBER 2019)
### REPLICATION FILE OF THE EMPIRICAL ANALYSIS
### JAN P. VOGLER
###########################

###########################
### INSTALL REQUIRED PACKAGES
###########################

install.packages("ggplot2")
install.packages("gridExtra")
install.packages("rddtools")
install.packages("MatchIt")
install.packages("rgenoud")
install.packages("BDSA")



###########################
### MAIN EMPIRICAL ANALYSIS IN THE ARTICLE
###########################

###########################
### READ DATASET (ALL PARTITIONS) AND DISPLAY SUMMARY STATISTICS
###########################

pol_00_complete=read.csv("pol_00_complete.csv", stringsAsFactors=F, encoding="UTF-8")
summary(pol_00_complete)



###########################
### SIMPLE DUMMY VARIABLE ANALYSIS (ALL PARTITIONS)
###########################

### EMPL. / POP. (LOG.)

lm0_dummy_1=lm(ln_empl_po ~ AU + RU + GER1918, data=pol_00_complete)
summary(lm0_dummy_1)

lm0_dummy_2=lm(ln_empl_po ~ AU + RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + gmina_wiej + gmina_mw, data=pol_00_complete)
summary(lm0_dummy_2)



### APP. / JOB (LOG.)

lm0_dummy_3=lm(ln_app_job ~ AU + RU + GER1918, data=pol_00_complete)
summary(lm0_dummy_3)

lm0_dummy_4=lm(ln_app_job ~ AU + RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop, data=pol_00_complete)
summary(lm0_dummy_4)



### ADVERT.

lm0_dummy_5_qp=glm(advert ~ AU + RU + GER1918, family="quasipoisson", data=pol_00_complete)
summary(lm0_dummy_5_qp)

lm0_dummy_6_qp=glm(advert ~ AU + RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop, family="quasipoisson", data=pol_00_complete)
summary(lm0_dummy_6_qp)



###########################
### COEFFICIENT PLOT: SIMPLE DUMMY VARIABLES (ALL PARTITIONS)
###########################

library(ggplot2)

### EMPL. / POP. (LOG.)

coefs = as.data.frame(summary(lm0_dummy_1)$coefficients[-c(1,5:11),1:2])
names(coefs)[2] = "se" 
coefs$vars = rownames(coefs)

p1 = ggplot(coefs, aes(vars, Estimate)) +
  ggtitle("EMPL. / POP.") +
  geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
  geom_errorbar(aes(ymin=Estimate - 1.96*se, ymax=Estimate + 1.96*se), 
                lwd=1.5, colour="firebrick", width=0) +
  geom_errorbar(aes(ymin=Estimate - 1.65*se, ymax=Estimate + 1.65*se), 
                lwd=2.5, colour="firebrick1", width=0) +
  geom_point(size=6, pch=21, fill="black") +
  theme_bw() + theme(text = element_text(size=20),
                     axis.text.x = element_text(angle=0),
                     plot.title = element_text(lineheight=.8, face="bold", hjust = 0.5))  +
  labs(x="Variables")



### EMPL. / POP. (LOG.) (WITH CONTROLS)

coefs = as.data.frame(summary(lm0_dummy_2)$coefficients[-c(1,5:12),1:2])
names(coefs)[2] = "se" 
coefs$vars = rownames(coefs)

p2 = ggplot(coefs, aes(vars, Estimate)) +
  ggtitle("With Controls") +
  geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
  geom_errorbar(aes(ymin=Estimate - 1.96*se, ymax=Estimate + 1.96*se), 
                lwd=1.5, colour="lightblue", width=0) +
  geom_errorbar(aes(ymin=Estimate - 1.65*se, ymax=Estimate + 1.65*se), 
                lwd=2.5, colour="blue", width=0) +
  geom_point(size=6, pch=21, fill="black") +
  theme_bw() + theme(text = element_text(size=20),
                     axis.text.x = element_text(angle=0),
                     plot.title = element_text(lineheight=.8, face="bold", hjust = 0.5))  +
  labs(x="Variables")



### APP. / JOB (LOG.)

coefs = as.data.frame(summary(lm0_dummy_3)$coefficients[-c(1,5:11),1:2])
names(coefs)[2] = "se" 
coefs$vars = rownames(coefs)

p3 = ggplot(coefs, aes(vars, Estimate)) +
  ggtitle("APP. / JOB") +
  geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
  geom_errorbar(aes(ymin=Estimate - 1.96*se, ymax=Estimate + 1.96*se), 
                lwd=1.5, colour="firebrick", width=0) +
  geom_errorbar(aes(ymin=Estimate - 1.65*se, ymax=Estimate + 1.65*se), 
                lwd=2.5, colour="firebrick1", width=0) +
  geom_point(size=6, pch=21, fill="black") +
  theme_bw() + theme(text = element_text(size=20),
                     axis.text.x = element_text(angle=0),
                     plot.title = element_text(lineheight=.8, face="bold", hjust = 0.5))  +
  labs(x="Variables")



### APP. / JOB (LOG.) (WITH CONTROLS)

coefs = as.data.frame(summary(lm0_dummy_4)$coefficients[-c(1,5:11),1:2])
names(coefs)[2] = "se" 
coefs$vars = rownames(coefs)

p4 = ggplot(coefs, aes(vars, Estimate)) +
  ggtitle("With Controls") +
  geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
  geom_errorbar(aes(ymin=Estimate - 1.96*se, ymax=Estimate + 1.96*se), 
                lwd=1.5, colour="lightblue", width=0) +
  geom_errorbar(aes(ymin=Estimate - 1.65*se, ymax=Estimate + 1.65*se), 
                lwd=2.5, colour="blue", width=0) +
  geom_point(size=6, pch=21, fill="black") +
  theme_bw() + theme(text = element_text(size=20),
                     axis.text.x = element_text(angle=0),
                     plot.title = element_text(lineheight=.8, face="bold", hjust = 0.5))  +
  labs(x="Variables")



###########################
### COMBINE ALL FOUR GRAPHS IN ONE PLOT
###########################

library(gridExtra)
grid.arrange(p1, p2, p3, p4, ncol=2)



###########################
### GEOGRAPHIC RD ANALYSIS
###########################

###########################
### PRUSSIA / RUSSIA COMPARISON
###########################

###########################
### READ DATASET (PRUSSIA / RUSSIA) AND DISPLAY SUMMARY STATISTICS
###########################

pol_01_pr_ru=read.csv("pol_01_pr_ru.csv", stringsAsFactors=F, encoding="UTF-8")
summary(pol_01_pr_ru)



###########################
### FULL SAMPLE
###########################

### EMPL. / POP. (LOG.) - AIR DISTANCE

lm1_distance_lnempl1=lm(ln_empl_po ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru)
summary(lm1_distance_lnempl1)

lm1_distance_lnempl2=lm(ln_empl_po ~ RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + gmina_wiej + gmina_mw + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru)
summary(lm1_distance_lnempl2)



### EMPL. / POP. (LOG.) - LAT. / LONG.

lm1_latlong_lnempl1=lm(ln_empl_po ~ RU + GER1918 + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru)
summary(lm1_latlong_lnempl1)

lm1_latlong_lnempl2=lm(ln_empl_po ~ RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + gmina_wiej + gmina_mw + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru)
summary(lm1_latlong_lnempl2)



### APP. / JOB (LOG.) - AIR DISTANCE

lm1_distance_lnapp1=lm(ln_app_job ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru)
summary(lm1_distance_lnapp1)

lm1_distance_lnapp2=lm(ln_app_job ~ RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru)
summary(lm1_distance_lnapp2)



### APP. / JOB (LOG.) - LAT. / LONG.

lm1_latlong_lnapp1=lm(ln_app_job ~ RU + GER1918 + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru)
summary(lm1_latlong_lnapp1)

lm1_latlong_lnapp2=lm(ln_app_job ~ RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru)
summary(lm1_latlong_lnapp2)



###########################
### FIND OPTIMAL BANDWIDTH
###########################

###########################
### CREATE SUBSETS WITH COMPLETE OBSERVATIONS OF DIFFERENT DEPENDENT VARIABLES
###########################

rdd1_data1 = pol_01_pr_ru[with(pol_01_pr_ru, complete.cases(ln_empl_po)), ]
rdd1_data2 = pol_01_pr_ru[with(pol_01_pr_ru, complete.cases(ln_app_job)), ]
rdd1_data3 = pol_01_pr_ru[with(pol_01_pr_ru, complete.cases(advert)), ]



###########################
### DV1: EMPL. / POP. (LOG.)
###########################

library(rddtools)

rdd1_obj1 = rdd_data(y = rdd1_data1$ln_empl_po, x = rdd1_data1$dist_pr_ru, covar = rdd1_data1, cutpoint = 0)
bw_ik_rdd1_obj1 = rdd_bw_ik(rdd1_obj1)
bw_ik_rdd1_obj1



###########################
### DV2: APP. / JOB (LOG.)
###########################

rdd1_obj2 = rdd_data(y = rdd1_data2$ln_app_job, x = rdd1_data2$dist_pr_ru, covar = rdd1_data2, cutpoint = 0)
bw_ik_rdd1_obj2 = rdd_bw_ik(rdd1_obj2)
bw_ik_rdd1_obj2



###########################
### DV3: ADVERT.
###########################

rdd1_obj3 = rdd_data(y = rdd1_data3$advert, x = rdd1_data3$dist_pr_ru, covar = rdd1_data3, cutpoint = 0)
bw_ik_rdd1_obj3 = rdd_bw_ik(rdd1_obj3)
bw_ik_rdd1_obj3



###########################
### CREATE BORDER SAMPLES
###########################

pol_01_pr_ru_bw100=pol_01_pr_ru[which(abs(pol_01_pr_ru$dist_pr_ru)<100),]
pol_01_pr_ru_bw125=pol_01_pr_ru[which(abs(pol_01_pr_ru$dist_pr_ru)<125),]
pol_01_pr_ru_bw135=pol_01_pr_ru[which(abs(pol_01_pr_ru$dist_pr_ru)<135),]
pol_01_pr_ru_bw150=pol_01_pr_ru[which(abs(pol_01_pr_ru$dist_pr_ru)<150),]
pol_01_pr_ru_bw155=pol_01_pr_ru[which(abs(pol_01_pr_ru$dist_pr_ru)<155),]
pol_01_pr_ru_bw175=pol_01_pr_ru[which(abs(pol_01_pr_ru$dist_pr_ru)<175),]
pol_01_pr_ru_bw200=pol_01_pr_ru[which(abs(pol_01_pr_ru$dist_pr_ru)<200),]
pol_01_pr_ru_bw225=pol_01_pr_ru[which(abs(pol_01_pr_ru$dist_pr_ru)<225),]
pol_01_pr_ru_bw250=pol_01_pr_ru[which(abs(pol_01_pr_ru$dist_pr_ru)<250),]
pol_01_pr_ru_bw260=pol_01_pr_ru[which(abs(pol_01_pr_ru$dist_pr_ru)<260),]



###########################
### BORDER SAMPLE REGRESSIONS
###########################

### EMPL. / POP. (LOG.)

lm1_rdd_bw100_lnempl=lm(ln_empl_po ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru_bw100)
summary(lm1_rdd_bw100_lnempl)

lm1_rdd_bw125_lnempl=lm(ln_empl_po ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru_bw125)
summary(lm1_rdd_bw125_lnempl)

lm1_rdd_bw135_lnempl=lm(ln_empl_po ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru_bw135)
summary(lm1_rdd_bw135_lnempl)

lm1_rdd_bw150_lnempl=lm(ln_empl_po ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru_bw150)
summary(lm1_rdd_bw150_lnempl)

lm1_rdd_bw175_lnempl=lm(ln_empl_po ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru_bw175)
summary(lm1_rdd_bw175_lnempl)

lm1_rdd_bw200_lnempl=lm(ln_empl_po ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru_bw200)
summary(lm1_rdd_bw200_lnempl)



### APP. / JOB (LOG.)

lm1_rdd_bw100_lnapp=lm(ln_app_job ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru_bw100)
summary(lm1_rdd_bw100_lnapp)

lm1_rdd_bw125_lnapp=lm(ln_app_job ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru_bw125)
summary(lm1_rdd_bw125_lnapp)

lm1_rdd_bw150_lnapp=lm(ln_app_job ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru_bw150)
summary(lm1_rdd_bw150_lnapp)

lm1_rdd_bw155_lnapp=lm(ln_app_job ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru_bw155)
summary(lm1_rdd_bw155_lnapp)

lm1_rdd_bw175_lnapp=lm(ln_app_job ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru_bw175)
summary(lm1_rdd_bw175_lnapp)

lm1_rdd_bw200_lnapp=lm(ln_app_job ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru_bw200)
summary(lm1_rdd_bw200_lnapp)



### ADVERT.

lm1_rdd_bw100_advt_qp=glm(advert ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, family = "quasipoisson", data=pol_01_pr_ru_bw100)
summary(lm1_rdd_bw100_advt_qp)

lm1_rdd_bw125_advt_qp=glm(advert ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, family = "quasipoisson", data=pol_01_pr_ru_bw125)
summary(lm1_rdd_bw125_advt_qp)

lm1_rdd_bw150_advt_qp=glm(advert ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, family = "quasipoisson", data=pol_01_pr_ru_bw150)
summary(lm1_rdd_bw150_advt_qp)

lm1_rdd_bw175_advt_qp=glm(advert ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, family = "quasipoisson", data=pol_01_pr_ru_bw175)
summary(lm1_rdd_bw175_advt_qp)

lm1_rdd_bw200_advt_qp=glm(advert ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, family = "quasipoisson", data=pol_01_pr_ru_bw200)
summary(lm1_rdd_bw200_advt_qp)

lm1_rdd_bw260_advt_qp=glm(advert ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, family = "quasipoisson", data=pol_01_pr_ru_bw260)
summary(lm1_rdd_bw260_advt_qp)



###########################
### AUSTRIA / RUSSIA COMPARISON
###########################

###########################
### READ DATASET (AUSTRIA / RUSSIA) AND DISPLAY SUMMARY STATISTICS
###########################

pol_02_au_ru=read.csv("pol_02_au_ru.csv", stringsAsFactors=F, encoding="UTF-8")
summary(pol_02_au_ru)



###########################
### FULL SAMPLE
###########################

### EMPL. / POP. (LOG.) - AIR DISTANCE

lm2_distance_lnempl1=lm(ln_empl_po ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru)
summary(lm2_distance_lnempl1)

lm2_distance_lnempl2=lm(ln_empl_po ~ RU + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + gmina_wiej + gmina_mw + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru)
summary(lm2_distance_lnempl2)



### EMPL. / POP. (LOG.) - LAT. / LONG.

lm2_latlong_lnempl1=lm(ln_empl_po ~ RU + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru)
summary(lm2_latlong_lnempl1)

lm2_latlong_lnempl2=lm(ln_empl_po ~ RU + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + gmina_wiej + gmina_mw + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru)
summary(lm2_latlong_lnempl2)



### APP. / JOB (LOG.) - AIR DISTANCE

lm2_distance_lnapp1=lm(ln_app_job ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru)
summary(lm2_distance_lnapp1)

lm2_distance_lnapp2=lm(ln_app_job ~ RU + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru)
summary(lm2_distance_lnapp2)



### APP. / JOB (LOG.) - LAT. / LONG.

lm2_latlong_lnapp1=lm(ln_app_job ~ RU + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru)
summary(lm2_latlong_lnapp1)

lm2_latlong_lnapp2=lm(ln_app_job ~ RU + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru)
summary(lm2_latlong_lnapp2)




###########################
### FIND OPTIMAL BANDWIDTH
###########################

###########################
### CREATE SUBSETS WITH COMPLETE OBSERVATIONS OF DIFFERENT DEPENDENT VARIABLES
###########################

rdd2_data1 = pol_02_au_ru[with(pol_02_au_ru, complete.cases(ln_empl_po)), ]
rdd2_data2 = pol_02_au_ru[with(pol_02_au_ru, complete.cases(ln_app_job)), ]
rdd2_data3 = pol_02_au_ru[with(pol_02_au_ru, complete.cases(advert)), ]



###########################
### DV1: EMPL. / POP. (LOG.)
###########################

rdd2_obj1 = rdd_data(y = rdd2_data1$ln_empl_po, x = rdd2_data1$dist_au_ru, covar = rdd2_data1, cutpoint = 0)
bw_ik_rdd2_obj1 = rdd_bw_ik(rdd2_obj1)
bw_ik_rdd2_obj1



###########################
### DV2: APP. / JOB (LOG.)
###########################

rdd2_obj2 = rdd_data(y = rdd2_data2$ln_app_job, x = rdd2_data2$dist_au_ru, covar = rdd2_data2, cutpoint = 0)
bw_ik_rdd2_obj2 = rdd_bw_ik(rdd2_obj2)
bw_ik_rdd2_obj2



###########################
### DV3: ADVERT.
###########################

rdd2_obj3 = rdd_data(y = rdd2_data3$advert, x = rdd2_data3$dist_au_ru, covar = rdd2_data3, cutpoint = 0)
bw_ik_rdd2_obj3 = rdd_bw_ik(rdd2_obj3)
bw_ik_rdd2_obj3



###########################
### CREATE BORDER SAMPLES
###########################

pol_02_au_ru_bw50=pol_02_au_ru[which(abs(pol_02_au_ru$dist_au_ru)<50),]
pol_02_au_ru_bw65=pol_02_au_ru[which(abs(pol_02_au_ru$dist_au_ru)<65),]
pol_02_au_ru_bw75=pol_02_au_ru[which(abs(pol_02_au_ru$dist_au_ru)<75),]
pol_02_au_ru_bw100=pol_02_au_ru[which(abs(pol_02_au_ru$dist_au_ru)<100),]
pol_02_au_ru_bw110=pol_02_au_ru[which(abs(pol_02_au_ru$dist_au_ru)<110),]
pol_02_au_ru_bw125=pol_02_au_ru[which(abs(pol_02_au_ru$dist_au_ru)<125),]
pol_02_au_ru_bw150=pol_02_au_ru[which(abs(pol_02_au_ru$dist_au_ru)<150),]
pol_02_au_ru_bw170=pol_02_au_ru[which(abs(pol_02_au_ru$dist_au_ru)<170),]
pol_02_au_ru_bw175=pol_02_au_ru[which(abs(pol_02_au_ru$dist_au_ru)<175),]



###########################
### BORDER SAMPLE REGRESSIONS
###########################

### EMPL. / POP. (LOG.)

lm2_rdd_bw50_lnempl=lm(ln_empl_po ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru_bw50)
summary(lm2_rdd_bw50_lnempl)

lm2_rdd_bw65_lnempl=lm(ln_empl_po ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru_bw65)
summary(lm2_rdd_bw65_lnempl)

lm2_rdd_bw75_lnempl=lm(ln_empl_po ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru_bw75)
summary(lm2_rdd_bw75_lnempl)

lm2_rdd_bw100_lnempl=lm(ln_empl_po ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru_bw100)
summary(lm2_rdd_bw100_lnempl)

lm2_rdd_bw125_lnempl=lm(ln_empl_po ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru_bw125)
summary(lm2_rdd_bw125_lnempl)

lm2_rdd_bw150_lnempl=lm(ln_empl_po ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru_bw150)
summary(lm2_rdd_bw150_lnempl)



### APP. / JOB (LOG.)

lm2_rdd_bw75_lnapp=lm(ln_app_job ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru_bw75)
summary(lm2_rdd_bw75_lnapp)

lm2_rdd_bw100_lnapp=lm(ln_app_job ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru_bw100)
summary(lm2_rdd_bw100_lnapp)

lm2_rdd_bw125_lnapp=lm(ln_app_job ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru_bw125)
summary(lm2_rdd_bw125_lnapp)

lm2_rdd_bw150_lnapp=lm(ln_app_job ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru_bw150)
summary(lm2_rdd_bw150_lnapp)

lm2_rdd_bw170_lnapp=lm(ln_app_job ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru_bw170)
summary(lm2_rdd_bw170_lnapp)

lm2_rdd_bw175_lnapp=lm(ln_app_job ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru_bw175)
summary(lm2_rdd_bw175_lnapp)



### ADVERT.

lm2_rdd_bw50_advt_qp=glm(advert ~ RU + dist_au_ru + RU*dist_au_ru, family = "quasipoisson", data=pol_02_au_ru_bw50)
summary(lm2_rdd_bw50_advt_qp)

lm2_rdd_bw75_advt_qp=glm(advert ~ RU + dist_au_ru + RU*dist_au_ru, family = "quasipoisson", data=pol_02_au_ru_bw75)
summary(lm2_rdd_bw75_advt_qp)

lm2_rdd_bw100_advt_qp=glm(advert ~ RU + dist_au_ru + RU*dist_au_ru, family = "quasipoisson", data=pol_02_au_ru_bw100)
summary(lm2_rdd_bw100_advt_qp)

lm2_rdd_bw110_advt_qp=glm(advert ~ RU + dist_au_ru + RU*dist_au_ru, family = "quasipoisson", data=pol_02_au_ru_bw110)
summary(lm2_rdd_bw110_advt_qp)

lm2_rdd_bw125_advt_qp=glm(advert ~ RU + dist_au_ru + RU*dist_au_ru, family = "quasipoisson", data=pol_02_au_ru_bw125)
summary(lm2_rdd_bw125_advt_qp)

lm2_rdd_bw150_advt_qp=glm(advert ~ RU + dist_au_ru + RU*dist_au_ru, family = "quasipoisson", data=pol_02_au_ru_bw150)
summary(lm2_rdd_bw150_advt_qp)



###########################
### PRUSSIA / AUSTRIA COMPARISON
###########################

###########################
### READ DATASET AND SHOW SUMMARY STATISTICS
###########################

pol_03_pr_au=read.csv("pol_03_pr_au.csv", stringsAsFactors=F, encoding="UTF-8")
summary(pol_03_pr_au)



###########################
### FULL SAMPLE
###########################

### EMPL. / POP. (LOG.) - AIR DISTANCE

lm3_distance_lnempl1=lm(ln_empl_po ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au)
summary(lm3_distance_lnempl1)

lm3_distance_lnempl2=lm(ln_empl_po ~ AU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + gmina_wiej + gmina_mw + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au)
summary(lm3_distance_lnempl2)



### EMPL. / POP. (LOG.) - LAT. / LONG.

lm3_latlong_lnempl1=lm(ln_empl_po ~ AU + GER1918 + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au)
summary(lm3_latlong_lnempl1)

lm3_latlong_lnempl2=lm(ln_empl_po ~ AU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + gmina_wiej + gmina_mw + academ_app + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au)
summary(lm3_latlong_lnempl2)



### APP. / JOB (LOG.) - AIR DISTANCE

lm3_distance_lnapp2=lm(ln_app_job ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au)
summary(lm3_distance_lnapp2)

lm3_distance_lnapp3=lm(ln_app_job ~ AU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au)
summary(lm3_distance_lnapp3)



### APP. / JOB (LOG.) - LAT. / LONG.

lm3_latlong_lnapp1=lm(ln_app_job ~ AU + GER1918 + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au)
summary(lm3_latlong_lnapp1)

lm3_latlong_lnapp2=lm(ln_app_job ~ AU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au)
summary(lm3_latlong_lnapp2)



###########################
### FIND OPTIMAL BANDWIDTH
###########################

###########################
### CREATE SUBSETS WITH COMPLETE OBSERVATIONS OF DIFFERENT DEPENDENT VARIABLES
###########################

rdd3_data1 = pol_03_pr_au[with(pol_03_pr_au, complete.cases(ln_empl_po)), ]
rdd3_data2 = pol_03_pr_au[with(pol_03_pr_au, complete.cases(ln_app_job)), ]
rdd3_data3 = pol_03_pr_au[with(pol_03_pr_au, complete.cases(advert)), ]



###########################
### DV1: EMPL. / POP. (LOG.)
###########################

rdd3_obj1 = rdd_data(y = rdd3_data1$ln_empl_po, x = rdd3_data1$dist_pr_au, covar = rdd3_data1, cutpoint = 0)
bw_ik_rdd3_obj1 = rdd_bw_ik(rdd3_obj1)
bw_ik_rdd3_obj1



###########################
### DV2: APP. / JOB (LOG.)
###########################

rdd3_obj2 = rdd_data(y = rdd3_data2$ln_app_job, x = rdd3_data2$dist_pr_au, covar = rdd3_data2, cutpoint = 0)
bw_ik_rdd3_obj2 = rdd_bw_ik(rdd3_obj2)
bw_ik_rdd3_obj2



###########################
### DV3: ADVERT.
###########################

rdd3_obj3 = rdd_data(y = rdd3_data3$advert, x = rdd3_data3$dist_pr_au, covar = rdd3_data3, cutpoint = 0)
bw_ik_rdd3_obj3 = rdd_bw_ik(rdd3_obj3)
bw_ik_rdd3_obj3



###########################
### CREATE BORDER SAMPLES
###########################

pol_03_pr_au_bw100=pol_03_pr_au[which(abs(pol_03_pr_au$dist_pr_au)<100),]
pol_03_pr_au_bw125=pol_03_pr_au[which(abs(pol_03_pr_au$dist_pr_au)<125),]
pol_03_pr_au_bw150=pol_03_pr_au[which(abs(pol_03_pr_au$dist_pr_au)<150),]
pol_03_pr_au_bw175=pol_03_pr_au[which(abs(pol_03_pr_au$dist_pr_au)<175),]
pol_03_pr_au_bw200=pol_03_pr_au[which(abs(pol_03_pr_au$dist_pr_au)<200),]
pol_03_pr_au_bw225=pol_03_pr_au[which(abs(pol_03_pr_au$dist_pr_au)<225),]
pol_03_pr_au_bw250=pol_03_pr_au[which(abs(pol_03_pr_au$dist_pr_au)<250),]
pol_03_pr_au_bw257=pol_03_pr_au[which(abs(pol_03_pr_au$dist_pr_au)<257),]
pol_03_pr_au_bw291=pol_03_pr_au[which(abs(pol_03_pr_au$dist_pr_au)<291),]



###########################
### BORDER SAMPLE REGRESSIONS
###########################

### EMPL. / POP. (LOG.) - AIR DISTANCE

lm3_rdd_bw100_lnempl=lm(ln_empl_po ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au_bw100)
summary(lm3_rdd_bw100_lnempl)

lm3_rdd_bw125_lnempl=lm(ln_empl_po ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au_bw125)
summary(lm3_rdd_bw125_lnempl)

lm3_rdd_bw150_lnempl=lm(ln_empl_po ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au_bw150)
summary(lm3_rdd_bw150_lnempl)

lm3_rdd_bw175_lnempl=lm(ln_empl_po ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au_bw175)
summary(lm3_rdd_bw175_lnempl)

lm3_rdd_bw200_lnempl=lm(ln_empl_po ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au_bw200)
summary(lm3_rdd_bw200_lnempl)



###########################
### MATCHING
###########################

###########################
### PRUSSIA / RUSSIA COMPARISON
###########################

library(MatchIt)



###########################
### DV1: EMPL. / POP. (LOG.)
###########################

pol_01_pr_ru_matching_dv1=pol_01_pr_ru[with(pol_01_pr_ru, complete.cases("ln_empl_po","RU","GER1918","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg"))]
pol_01_pr_ru_matching_dv1=pol_01_pr_ru_matching_dv1[,c("ln_empl_po","RU","GER1918","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg")]
pol_01_pr_ru_matching_dv1=na.omit(pol_01_pr_ru_matching_dv1)
set.seed(06112019) # JUNE 11, 2019
matching1_dv1 = matchit(RU ~ ln_revenue + ln_pop_den + migr_avg + unempl_avg + city_powia + GER1918, data=pol_01_pr_ru_matching_dv1, method = "genetic", ratio = 1)
summary(matching1_dv1)
matching1.data.dv1 = match.data(matching1_dv1)
lm_matching_data1_dv1 = lm(ln_empl_po ~ RU, data = matching1.data.dv1)
summary(lm_matching_data1_dv1)



###########################
### DV2: APP. / JOB (LOG.)
###########################

pol_01_pr_ru_matching_dv2=pol_01_pr_ru[with(pol_01_pr_ru, complete.cases("ln_app_job","RU","GER1918","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg","lnpop"))]
pol_01_pr_ru_matching_dv2=pol_01_pr_ru_matching_dv2[,c("ln_app_job","RU","GER1918","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg","lnpop")]
pol_01_pr_ru_matching_dv2=na.omit(pol_01_pr_ru_matching_dv2)
set.seed(06112019) # JUNE 11, 2019
matching1_dv2 = matchit(RU ~ ln_revenue + ln_pop_den + migr_avg + unempl_avg + city_powia + GER1918 + lnpop, data=pol_01_pr_ru_matching_dv2, method = "genetic", ratio = 1)
summary(matching1_dv2)
matching1.data.dv2 = match.data(matching1_dv2)
lm_matching_data1_dv2 = lm(ln_app_job ~ RU, data = matching1.data.dv2)
summary(lm_matching_data1_dv2)



###########################
### DV3: ADVERT.
###########################

pol_01_pr_ru_matching_dv3=pol_01_pr_ru[with(pol_01_pr_ru, complete.cases("advert","RU","GER1918","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg","lnpop"))]
pol_01_pr_ru_matching_dv3=pol_01_pr_ru_matching_dv3[,c("advert","RU","GER1918","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg","lnpop")]
pol_01_pr_ru_matching_dv3=na.omit(pol_01_pr_ru_matching_dv3)
set.seed(06112019) # JUNE 11, 2019
matching1_dv3 = matchit(RU ~ ln_revenue + ln_pop_den + migr_avg + unempl_avg + city_powia + GER1918 + lnpop, data=pol_01_pr_ru_matching_dv3, method = "genetic", ratio = 1)
summary(matching1_dv3)
matching1.data.dv3 = match.data(matching1_dv3)
lm_matching_data1_dv3_qp = glm(advert ~ RU, family = "quasipoisson", data = matching1.data.dv3)
summary(lm_matching_data1_dv3_qp)



###########################
### AUSTRIA / RUSSIA COMPARISON
###########################

###########################
### DV1: EMPL. / POP. (LOG.)
###########################

pol_02_au_ru_matching_dv1=pol_02_au_ru[with(pol_02_au_ru, complete.cases("ln_empl_po","RU","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg"))]
pol_02_au_ru_matching_dv1=pol_02_au_ru_matching_dv1[,c("ln_empl_po","RU","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg")]
pol_02_au_ru_matching_dv1=na.omit(pol_02_au_ru_matching_dv1)
set.seed(06112019) # JUNE 11, 2019
matching2_dv1 = matchit(RU ~ ln_revenue + ln_pop_den + migr_avg + unempl_avg + city_powia, data=pol_02_au_ru_matching_dv1, method = "genetic", ratio = 2)
summary(matching2_dv1)
matching2.data.dv1 = match.data(matching2_dv1)
lm_matching_data2_dv1 = lm(ln_empl_po ~ RU, data = matching2.data.dv1)
summary(lm_matching_data2_dv1)



###########################
### DV2: APP. / JOB (LOG.)
###########################

pol_02_au_ru_matching_dv2=pol_02_au_ru[with(pol_02_au_ru, complete.cases("ln_app_job","RU","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg","lnpop"))]
pol_02_au_ru_matching_dv2=pol_02_au_ru_matching_dv2[,c("ln_app_job","RU","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg","lnpop")]
pol_02_au_ru_matching_dv2=na.omit(pol_02_au_ru_matching_dv2)
set.seed(06112019) # JUNE 11, 2019
matching2_dv2 = matchit(RU ~ ln_revenue + ln_pop_den + migr_avg + unempl_avg + city_powia + lnpop, data=pol_02_au_ru_matching_dv2, method = "genetic", ratio = 1)
summary(matching2_dv2)
matching2.data.dv2 = match.data(matching2_dv2)
lm_matching_data2_dv2 = lm(ln_app_job ~ RU, data = matching2.data.dv2)
summary(lm_matching_data2_dv2)



###########################
### DV3: ADVERT.
###########################

pol_02_au_ru_matching_dv3=pol_02_au_ru[with(pol_02_au_ru, complete.cases("advert","RU","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg","lnpop"))]
pol_02_au_ru_matching_dv3=pol_02_au_ru_matching_dv3[,c("advert","RU","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg","lnpop")]
pol_02_au_ru_matching_dv3=na.omit(pol_02_au_ru_matching_dv3)
set.seed(06112019) # JUNE 11, 2019
matching2_dv3 = matchit(RU ~ ln_revenue + ln_pop_den + migr_avg + unempl_avg + city_powia + lnpop, data=pol_02_au_ru_matching_dv3, method = "genetic", ratio = 2)
summary(matching2_dv3)
matching2.data.dv3 = match.data(matching2_dv3)
lm_matching_data2_dv3_qp = glm(advert ~ RU, family = "quasipoisson", data = matching2.data.dv3)
summary(lm_matching_data2_dv3_qp)



###########################
### PRUSSIA / AUSTRIA COMPARISON
###########################

###########################
### DV1: EMPL. / POP. (LOG.)
###########################

pol_03_pr_au_matching_dv1=pol_03_pr_au[with(pol_03_pr_au, complete.cases("ln_empl_po","AU","GER1918","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg"))]
pol_03_pr_au_matching_dv1=pol_03_pr_au_matching_dv1[,c("ln_empl_po","AU","GER1918","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg")]
pol_03_pr_au_matching_dv1=na.omit(pol_03_pr_au_matching_dv1)
set.seed(06112019) # JUNE 11, 2019
matching3_dv1 = matchit(AU ~ ln_revenue + ln_pop_den + migr_avg + unempl_avg + city_powia + GER1918, data=pol_03_pr_au_matching_dv1, method = "genetic", ratio = 1)
summary(matching3_dv1)
matching3.data.dv1 = match.data(matching3_dv1)
lm_matching_data3_dv1 = lm(ln_empl_po ~ AU, data = matching3.data.dv1)
summary(lm_matching_data3_dv1)



###########################
### DV2: APP. / JOB (LOG.)
###########################

pol_03_pr_au_matching_dv2=pol_03_pr_au[with(pol_03_pr_au, complete.cases("ln_app_job","AU","GER1918","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg","lnpop"))]
pol_03_pr_au_matching_dv2=pol_03_pr_au_matching_dv2[,c("ln_app_job","AU","GER1918","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg","lnpop")]
summary(pol_03_pr_au_matching_dv2)
pol_03_pr_au_matching_dv2=na.omit(pol_03_pr_au_matching_dv2)
set.seed(06112019) # JUNE 11, 2019
matching3_dv2 = matchit(AU ~ ln_revenue + ln_pop_den + migr_avg + unempl_avg + city_powia + GER1918 + lnpop, data=pol_03_pr_au_matching_dv2, method = "genetic", ratio = 1)
summary(matching3_dv2)
matching3.data.dv2 = match.data(matching3_dv2)
lm_matching_data3_dv2 = lm(ln_app_job ~ AU, data = matching3.data.dv2)
summary(lm_matching_data3_dv2)



###########################
### DV3: ADVERT.
###########################

pol_03_pr_au_matching_dv3=pol_03_pr_au[with(pol_03_pr_au, complete.cases("advert","AU","GER1918","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg","lnpop"))]
pol_03_pr_au_matching_dv3=pol_03_pr_au_matching_dv3[,c("advert","AU","GER1918","ln_revenue","ln_pop_den","city_powia","migr_avg","unempl_avg","lnpop")]
pol_03_pr_au_matching_dv3=na.omit(pol_03_pr_au_matching_dv3)
set.seed(06112019) # JUNE 11, 2019
matching3_dv3 = matchit(AU ~ ln_revenue + ln_pop_den + migr_avg + unempl_avg + city_powia + GER1918 + lnpop, data=pol_03_pr_au_matching_dv3, method = "genetic", ratio = 1)
summary(matching3_dv3)
matching3.data.dv3 = match.data(matching3_dv3)
lm_matching_data3_dv3_qp = glm(advert ~ AU, family = "quasipoisson", data = matching3.data.dv3)
summary(lm_matching_data3_dv3_qp)



###########################
### ADDITIONAL EMPIRICAL ANALYSIS IN THE SUPPLEMENTARY MATERIAL
###########################

###########################
### PRE-TREATMENT CHARACTERISTIC COMPARISON
###########################

pol_04_pre_tr=read.csv("pol_04_pre_tr.csv", stringsAsFactors=F, encoding="UTF-8")
summary(pol_04_pre_tr)

library(BSDA)



###########################
### PRUSSIAN VS. NON-PRUSSIAN TOWNS
###########################

t.test(pol_04_pre_tr$pop1400[pol_04_pre_tr$PR==1], pol_04_pre_tr$pop1400[pol_04_pre_tr$PR==0])
z.test(pol_04_pre_tr$major_trad[pol_04_pre_tr$PR==1], pol_04_pre_tr$major_trad[pol_04_pre_tr$PR==0], sigma.x=sd(pol_04_pre_tr$major_trad[pol_04_pre_tr$PR==1]), sigma.y=sd(pol_04_pre_tr$major_trad[pol_04_pre_tr$PR==0]))
z.test(pol_04_pre_tr$see_of_ecc[pol_04_pre_tr$PR==1], pol_04_pre_tr$see_of_ecc[pol_04_pre_tr$PR==0], sigma.x=sd(pol_04_pre_tr$see_of_ecc[pol_04_pre_tr$PR==1]), sigma.y=sd(pol_04_pre_tr$see_of_ecc[pol_04_pre_tr$PR==0]))



###########################
### AUSTRIAN VS. NON-AUSTRIAN TOWNS
###########################

t.test(pol_04_pre_tr$pop1400[pol_04_pre_tr$AU==1], pol_04_pre_tr$pop1400[pol_04_pre_tr$AU==0])
z.test(pol_04_pre_tr$major_trad[pol_04_pre_tr$AU==1], pol_04_pre_tr$major_trad[pol_04_pre_tr$AU==0], sigma.x=sd(pol_04_pre_tr$major_trad[pol_04_pre_tr$AU==1]), sigma.y=sd(pol_04_pre_tr$major_trad[pol_04_pre_tr$AU==0]))
z.test(pol_04_pre_tr$see_of_ecc[pol_04_pre_tr$AU==1], pol_04_pre_tr$see_of_ecc[pol_04_pre_tr$AU==0], sigma.x=sd(pol_04_pre_tr$see_of_ecc[pol_04_pre_tr$AU==1]), sigma.y=sd(pol_04_pre_tr$see_of_ecc[pol_04_pre_tr$AU==0]))



###########################
### RUSSIAN VS. NON-RUSSIAN TOWNS
###########################

t.test(pol_04_pre_tr$pop1400[pol_04_pre_tr$RU==1], pol_04_pre_tr$pop1400[pol_04_pre_tr$RU==0])
z.test(pol_04_pre_tr$major_trad[pol_04_pre_tr$RU==1], pol_04_pre_tr$major_trad[pol_04_pre_tr$RU==0], sigma.x=sd(pol_04_pre_tr$major_trad[pol_04_pre_tr$RU==1]), sigma.y=sd(pol_04_pre_tr$major_trad[pol_04_pre_tr$RU==0]))
z.test(pol_04_pre_tr$see_of_ecc[pol_04_pre_tr$RU==1], pol_04_pre_tr$see_of_ecc[pol_04_pre_tr$RU==0], sigma.x=sd(pol_04_pre_tr$see_of_ecc[pol_04_pre_tr$RU==1]), sigma.y=sd(pol_04_pre_tr$see_of_ecc[pol_04_pre_tr$RU==0]))



###########################
### MERITOCRACY AND EFFICIENCY
###########################

lm0_eff_1=lm(ln_time_ve ~ ln_app_job, data=pol_00_complete)
summary(lm0_eff_1)

lm0_eff_2=lm(ln_time_ve ~ ln_app_job + lnvehicle + WARSAW_DIS + ln_revenue + ln_pop_den + migr_avg + unempl_avg, data=pol_00_complete)
summary(lm0_eff_2)

lm0_eff_3=lm(ln_empl_po ~ ln_app_job, data=pol_00_complete)
summary(lm0_eff_3)

lm0_eff_4=lm(ln_empl_po ~ ln_app_job + ln_revenue + ln_pop_den + migr_avg + unempl_avg + city_powia + gmina_wiej + gmina_mw, data=pol_00_complete)
summary(lm0_eff_4)



###########################
### RELATIVE EMPLOYEES AND SPEED
###########################

lm0_eff_5=lm(ln_time_ve ~ ln_empl_po, data=pol_00_complete)
summary(lm0_eff_5)

lm0_eff_6=lm(ln_time_ve ~ ln_empl_po + lnvehicle + WARSAW_DIS + ln_revenue + ln_pop_den + migr_avg + unempl_avg, data=pol_00_complete)
summary(lm0_eff_6)

lm0_eff_7=lm(ln_difftim ~ ln_empl_po, data=pol_00_complete)
summary(lm0_eff_7)

lm0_eff_8=lm(ln_difftim ~ ln_empl_po + ln_revenue + ln_pop_den + migr_avg + unempl_avg + lnpop + city_powia, data=pol_00_complete)
summary(lm0_eff_8)



###########################
### ADDITIONAL SIMPLE DUMMY VARIABLE ANALYSIS (AT OPTIMAL BANDWIDTHS)
###########################

###########################
### PRUSSIA / RUSSIA COMPARISON
###########################

pr_ru_limitedbw_emp=lm(ln_empl_po ~ RU + GER1918 + city_powia, data=pol_01_pr_ru_bw135)
summary(pr_ru_limitedbw_emp)

pr_ru_limitedbw_app=lm(ln_app_job ~ RU + GER1918 + city_powia, data=pol_01_pr_ru_bw155)
summary(pr_ru_limitedbw_app)

pr_ru_limitedbw_advert_qp=glm(advert ~ RU + GER1918 + city_powia, family = "quasipoisson", data=pol_01_pr_ru_bw200)
summary(pr_ru_limitedbw_advert_qp)



###########################
### AUSTRIA / RUSSIA COMPARISON
###########################

au_ru_limitedbw_emp=lm(ln_empl_po ~ RU + city_powia, data=pol_02_au_ru_bw65)
summary(au_ru_limitedbw_emp)

au_ru_limitedbw_app=lm(ln_app_job ~ RU + city_powia, data=pol_02_au_ru_bw170)
summary(au_ru_limitedbw_app)

au_ru_limitedbw_advert_qp=glm(advert ~ RU + city_powia, family = "quasipoisson", data=pol_02_au_ru_bw110)
summary(au_ru_limitedbw_advert_qp)



###########################
### PRUSSIA / AUSTRIA COMPARISON
###########################

pr_au_limitedbw_emp=lm(ln_empl_po ~ AU + GER1918 + city_powia, data=pol_03_pr_au_bw200)
summary(pr_au_limitedbw_emp)

pr_au_limitedbw_app=lm(ln_app_job ~ AU + GER1918 + city_powia, data=pol_03_pr_au_bw291)
summary(pr_au_limitedbw_app)

pr_au_limitedbw_advert_qp=glm(advert ~ AU + GER1918 + city_powia, family = "quasipoisson", data=pol_03_pr_au_bw257)
summary(pr_au_limitedbw_advert_qp)



###########################
### CORRECTING P-VALUES FOR MULTIPLE COMPARISONS
###########################

p_val_old_limitedbw_all=c(summary(pr_ru_limitedbw_emp)$coefficients[,4][2],summary(pr_ru_limitedbw_app)$coefficients[,4][2],summary(pr_ru_limitedbw_advert_qp)$coefficients[,4][2],summary(au_ru_limitedbw_emp)$coefficients[,4][2],summary(au_ru_limitedbw_app)$coefficients[,4][2],summary(au_ru_limitedbw_advert_qp)$coefficients[,4][2],summary(pr_au_limitedbw_emp)$coefficients[,4][2],summary(pr_au_limitedbw_app)$coefficients[,4][2],summary(pr_au_limitedbw_advert_qp)$coefficients[,4][2])
p_val_new_limitedbw_all=p.adjust(p_val_old_limitedbw_all, method = "holm", n = 9)
p_val_new_limitedbw_all # CORRECTED P-VALUES



###########################
### IMPERIAL LEGACIES IN OTHER DIMENSIONS
###########################

lm_pt_bias_rev=lm(ln_revenue ~ AU + RU + GER1918 + city_powia, data=pol_00_complete)
summary(lm_pt_bias_rev)

lm_pt_bias_pop_den=lm(ln_pop_den ~ AU + RU + GER1918 + city_powia, data=pol_00_complete)
summary(lm_pt_bias_pop_den)

lm_pt_bias_unempl_avg=lm(unempl_avg ~ AU + RU + GER1918 + city_powia, data=pol_00_complete)
summary(lm_pt_bias_unempl_avg)



###########################
### ADDITIONAL GEOGRAPHIC RD ANALYSES
###########################

###########################
### PRUSSIA / RUSSIA COMPARISON
###########################

###########################
### FULL SAMPLE
###########################

### ADVERT. - AIR DISTANCE

lm1_distance_advt1_qp=glm(advert ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, family = "quasipoisson", data=pol_01_pr_ru)
summary(lm1_distance_advt1_qp)

lm1_distance_advt2_qp=glm(advert ~ RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + dist_pr_ru + RU*dist_pr_ru, family = "quasipoisson", data=pol_01_pr_ru)
summary(lm1_distance_advt2_qp)



### ADVERT. - LAT. / LONG.

lm1_latlong_advt1_qp=glm(advert ~ RU + GER1918 + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_ru + RU*dist_pr_ru, family = "quasipoisson", data=pol_01_pr_ru)
summary(lm1_latlong_advt1_qp)

lm1_latlong_advt2_qp=glm(advert ~ RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_ru + RU*dist_pr_ru, family = "quasipoisson", data=pol_01_pr_ru)
summary(lm1_latlong_advt2_qp)



###########################
### RD GRAPHS
###########################

rdd1_obj1_lm = rdd_reg_lm(rdd_object = rdd1_obj1, covariates="GER1918", bw=200)
rdd1_obj2_lm = rdd_reg_lm(rdd_object = rdd1_obj2, covariates="GER1918", bw=200)
rdd1_obj3_lm = rdd_reg_lm(rdd_object = rdd1_obj3, covariates="GER1918", bw=200)



### DENSITY TESTS

dens_test(rdd1_obj1)
dens_test(rdd1_obj2)
dens_test(rdd1_obj3)



### SENSITIVITY TESTS

plotSensi(rdd1_obj1_lm, from = 100, to = 200, order=1:2, level = 0.9)
plotSensi(rdd1_obj2_lm, from = 100, to = 200, order=1:2, level = 0.9)
plotSensi(rdd1_obj3_lm, from = 100, to = 200, order=1:2, level = 0.9)



### PLACEBO TESTS

plotPlacebo(rdd1_obj1_lm, from=0.15, by=10)
plotPlacebo(rdd1_obj2_lm, from=0.15, by=10)
plotPlacebo(rdd1_obj3_lm, from=0.15, by=10)



###########################
### AUSTRIA / RUSSIA COMPARISON
###########################

###########################
### FULL SAMPLE
###########################

### ADVERT. - AIR DISTANCE

lm2_distance_advt1_qp=glm(advert ~ RU + dist_au_ru + RU*dist_au_ru, family = "quasipoisson", data=pol_02_au_ru)
summary(lm2_distance_advt1_qp)

lm2_distance_advt2_qp=glm(advert ~ RU + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + dist_au_ru + RU*dist_au_ru, family = "quasipoisson", data=pol_02_au_ru)
summary(lm2_distance_advt2_qp)



### ADVERT. - LAT. / LONG.

lm2_latlong_advt1_qp=glm(advert ~ RU + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_au_ru + RU*dist_au_ru, family = "quasipoisson", data=pol_02_au_ru)
summary(lm2_latlong_advt1_qp)

lm2_latlong_advt2_qp=glm(advert ~ RU + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_au_ru + RU*dist_au_ru, family = "quasipoisson", data=pol_02_au_ru)
summary(lm2_latlong_advt2_qp)



###########################
### RD GRAPHS
###########################

rdd2_obj1_lm = rdd_reg_lm(rdd_object = rdd2_obj1, bw=150)
rdd2_obj2_lm = rdd_reg_lm(rdd_object = rdd2_obj2, bw=175)
rdd2_obj3_lm = rdd_reg_lm(rdd_object = rdd2_obj3, bw=150)



### DENSITY TESTS

dens_test(rdd2_obj1)
dens_test(rdd2_obj2)
dens_test(rdd2_obj3)



### SENSITIVITY TESTS

plotSensi(rdd2_obj1_lm, from = 50, to = 150, order=1:2, level = 0.9)
plotSensi(rdd2_obj2_lm, from = 75, to = 175, order=1:2, level = 0.9)
plotSensi(rdd2_obj3_lm, from = 50, to = 150, order=1:2, level = 0.9)



### PLACEBO TESTS

plotPlacebo(rdd2_obj1_lm, from=0.15, by=10)
plotPlacebo(rdd2_obj2_lm, from=0.15, by=10)
plotPlacebo(rdd2_obj3_lm, from=0.15, by=10)



###########################
### PRUSSIA / AUSTRIA COMPARISON
###########################

###########################
### FULL SAMPLE
###########################

### ADVERT. - AIR DISTANCE

lm3_distance_advt1_qp=glm(advert ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, family = "quasipoisson", data=pol_03_pr_au)
summary(lm3_distance_advt1_qp)

lm3_distance_advt2_qp=glm(advert ~ AU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + dist_pr_au + AU*dist_pr_au, family = "quasipoisson", data=pol_03_pr_au)
summary(lm3_distance_advt2_qp)



### ADVERT. - LAT. / LONG.

lm3_latlong_advt1_qp=glm(advert ~ AU + GER1918 + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_au + AU*dist_pr_au, family = "quasipoisson", data=pol_03_pr_au)
summary(lm3_latlong_advt1_qp)

lm3_latlong_advt2_qp=glm(advert ~ AU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_au + AU*dist_pr_au, family = "quasipoisson", data=pol_03_pr_au)
summary(lm3_latlong_advt2_qp)



###########################
### BORDER SAMPLE REGRESSIONS
###########################

### APP. / JOB (LOG.)

lm3_rdd_bw100_lnapp=lm(ln_app_job ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au_bw100)
summary(lm3_rdd_bw100_lnapp)

lm3_rdd_bw125_lnapp=lm(ln_app_job ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au_bw125)
summary(lm3_rdd_bw125_lnapp)

lm3_rdd_bw150_lnapp=lm(ln_app_job ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au_bw150)
summary(lm3_rdd_bw150_lnapp)

lm3_rdd_bw175_lnapp=lm(ln_app_job ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au_bw175)
summary(lm3_rdd_bw175_lnapp)

lm3_rdd_bw200_lnapp=lm(ln_app_job ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au_bw200)
summary(lm3_rdd_bw200_lnapp)

lm3_rdd_bw291_lnapp=lm(ln_app_job ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au_bw291)
summary(lm3_rdd_bw291_lnapp)


### ADVERT.

lm3_rdd_bw100_advt_qp=glm(advert ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, family = "quasipoisson", data=pol_03_pr_au_bw100)
summary(lm3_rdd_bw100_advt_qp)

lm3_rdd_bw125_advt_qp=glm(advert ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, family = "quasipoisson", data=pol_03_pr_au_bw125)
summary(lm3_rdd_bw125_advt_qp)

lm3_rdd_bw150_advt_qp=glm(advert ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, family = "quasipoisson", data=pol_03_pr_au_bw150)
summary(lm3_rdd_bw150_advt_qp)

lm3_rdd_bw175_advt_qp=glm(advert ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, family = "quasipoisson", data=pol_03_pr_au_bw175)
summary(lm3_rdd_bw175_advt_qp)

lm3_rdd_bw200_advt_qp=glm(advert ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, family = "quasipoisson", data=pol_03_pr_au_bw200)
summary(lm3_rdd_bw200_advt_qp)

lm3_rdd_bw257_advt_qp=glm(advert ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, family = "quasipoisson", data=pol_03_pr_au_bw257)
summary(lm3_rdd_bw257_advt_qp)



###########################
### RD GRAPHS
###########################

rdd3_obj1_lm = rdd_reg_lm(rdd_object = rdd3_obj1, covariates="GER1918", bw=200)
rdd3_obj2_lm = rdd_reg_lm(rdd_object = rdd3_obj2, covariates="GER1918", bw=200)
rdd3_obj3_lm = rdd_reg_lm(rdd_object = rdd3_obj3, covariates="GER1918", bw=200)



### DENSITY TESTS

dens_test(rdd3_obj1)
dens_test(rdd3_obj2)
dens_test(rdd3_obj3)



### SENSITIVITY TESTS

plotSensi(rdd3_obj1_lm, from = 100, to = 200, order=1:2, level = 0.9)
plotSensi(rdd3_obj2_lm, from = 100, to = 200, order=1:2, level = 0.9)
plotSensi(rdd3_obj3_lm, from = 100, to = 200, order=1:2, level = 0.9)



### PLACEBO TESTS

plotPlacebo(rdd3_obj1_lm, from=0.15, by=10)
plotPlacebo(rdd3_obj2_lm, from=0.15, by=10)
plotPlacebo(rdd3_obj3_lm, from=0.15, by=10)



###########################################
### EXTENSIONS OF THE EMPIRICAL ANALYSIS
##########################################

###########################
### EXTENSION 1 - VARIATION WITHIN PRESENT-DAY VOIVODESHIPS
###########################

###########################
### PRUSSIA / RUSSIA COMPARISON
###########################

pol_01_pr_ru_subset=pol_01_pr_ru[which(pol_01_pr_ru$voivod=="kujawsko pomorskie" | pol_01_pr_ru$voivod=="wielkopolskie" | pol_01_pr_ru$voivod=="slaskie"),]

lm_new_ext_01=lm(ln_empl_po ~ RU + GER1918 + as.factor(voivod)-1, data=pol_01_pr_ru_subset)
summary(lm_new_ext_01)

lm_new_ext_02=lm(ln_app_job ~ RU + GER1918 + as.factor(voivod)-1, data=pol_01_pr_ru_subset)
summary(lm_new_ext_02)

lm_new_ext_03_qp=glm(advert ~ RU + GER1918 + as.factor(voivod)-1, family="quasipoisson", data=pol_01_pr_ru_subset)
summary(lm_new_ext_03_qp)



###########################
### AUSTRIA / RUSSIA COMPARISON
###########################

pol_02_au_ru_subset=pol_02_au_ru[which(pol_02_au_ru$voivod=="malopolskie" | pol_02_au_ru$voivod=="slaskie"),]

lm_new_ext_04=lm(ln_empl_po ~ RU + as.factor(voivod)-1, data=pol_02_au_ru_subset)
summary(lm_new_ext_04)

lm_new_ext_05=lm(ln_app_job ~ RU + as.factor(voivod)-1, data=pol_02_au_ru_subset)
summary(lm_new_ext_05)

lm_new_ext_06_qp=glm(advert ~ RU + as.factor(voivod)-1, family="quasipoisson", data=pol_02_au_ru_subset)
summary(lm_new_ext_06_qp)



###########################
### PRUSSIA / AUSTRIA COMPARISON
###########################

pol_03_pr_au_subset=pol_03_pr_au[which(pol_03_pr_au$voivod=="slaskie"),]

lm_new_ext_07=lm(ln_empl_po ~ AU, data=pol_03_pr_au_subset)
summary(lm_new_ext_07)

lm_new_ext_08=lm(ln_app_job ~ AU, data=pol_03_pr_au_subset)
summary(lm_new_ext_08)

lm_new_ext_09_qp=glm(advert ~ AU, family="quasipoisson", data=pol_03_pr_au_subset)
summary(lm_new_ext_09_qp)



###########################
### EXTENSION 2 - MAYORAL POLITICAL AFFILIATION AND REGIONAL GDP
###########################

### EMPL. / POP. (LOG.)

lm0_ext_dummy_1=lm(ln_empl_po ~ AU + RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + mayor_sld + mayor_po + mayor_psl + mayor_pis + ln_GDP + gmina_wiej + gmina_mw, data=pol_00_complete)
summary(lm0_ext_dummy_1)



### APP. / JOB (LOG.)

lm0_ext_dummy_2=lm(ln_app_job ~ AU + RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + mayor_sld + mayor_po + mayor_psl + mayor_pis + ln_GDP + lnpop, data=pol_00_complete)
summary(lm0_ext_dummy_2)



### ADVERT.

lm0_ext_dummy_3_qp=glm(advert ~ AU + RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + mayor_sld + mayor_po + mayor_psl + mayor_pis + ln_GDP + lnpop, family="quasipoisson", data=pol_00_complete)
summary(lm0_ext_dummy_3_qp)



###########################
### EXTENSION 3 - INCLUDING DISTANCE WEIGHTS
###########################

###########################
### PRUSSIA / RUSSIA COMPARISON
###########################

### EMPL. / POP. (LOG.) - AIR DISTANCE

lm1_dist_weights_lnempl1=lm(ln_empl_po ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru, weight=pol_01_pr_ru$dist_weights)
summary(lm1_dist_weights_lnempl1)

lm1_dist_weights_lnempl2=lm(ln_empl_po ~ RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + gmina_wiej + gmina_mw + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru, weight=pol_01_pr_ru$dist_weights)
summary(lm1_dist_weights_lnempl2)



### EMPL. / POP. (LOG.) - LAT. / LONG.

lm1_latlong_weights_lnempl1=lm(ln_empl_po ~ RU + GER1918 + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru, weight=pol_01_pr_ru$dist_weights)
summary(lm1_latlong_weights_lnempl1)

lm1_latlong_weights_lnempl2=lm(ln_empl_po ~ RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + gmina_wiej + gmina_mw + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru, weight=pol_01_pr_ru$dist_weights)
summary(lm1_latlong_weights_lnempl2)



### APP. / JOB (LOG.) - AIR DISTANCE

lm1_dist_weights_lnapp1=lm(ln_app_job ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru, weight=pol_01_pr_ru$dist_weights)
summary(lm1_dist_weights_lnapp1)

lm1_dist_weights_lnapp2=lm(ln_app_job ~ RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru, weight=pol_01_pr_ru$dist_weights)
summary(lm1_dist_weights_lnapp2)



### APP. / JOB (LOG.) - LAT. / LONG.

lm1_latlong_weights_lnapp1=lm(ln_app_job ~ RU + GER1918 + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru, weight=pol_01_pr_ru$dist_weights)
summary(lm1_latlong_weights_lnapp1)

lm1_latlong_weights_lnapp2=lm(ln_app_job ~ RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_ru + RU*dist_pr_ru, data=pol_01_pr_ru, weight=pol_01_pr_ru$dist_weights)
summary(lm1_latlong_weights_lnapp2)



### ADVERT. - AIR DISTANCE

lm1_dist_weights_advt1_qp=glm(advert ~ RU + GER1918 + dist_pr_ru + RU*dist_pr_ru, family="quasipoisson", data=pol_01_pr_ru, weight=pol_01_pr_ru$dist_weights)
summary(lm1_dist_weights_advt1_qp)

lm1_dist_weights_advt2_qp=glm(advert ~ RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + dist_pr_ru + RU*dist_pr_ru, family="quasipoisson", data=pol_01_pr_ru, weight=pol_01_pr_ru$dist_weights)
summary(lm1_dist_weights_advt2_qp)



### ADVERT. - LAT. / LONG.

lm1_latlong_weights_advt1_qp=glm(advert ~ RU + GER1918 + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_ru + RU*dist_pr_ru, family="quasipoisson", data=pol_01_pr_ru, weight=pol_01_pr_ru$dist_weights)
summary(lm1_latlong_weights_advt1_qp)

lm1_latlong_weights_advt2_qp=glm(advert ~ RU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_ru + RU*dist_pr_ru, family="quasipoisson", data=pol_01_pr_ru, weight=pol_01_pr_ru$dist_weights)
summary(lm1_latlong_weights_advt2_qp)



###########################
### AUSTRIA / RUSSIA COMPARISON
###########################

### EMPL. / POP. (LOG.) - AIR DISTANCE

lm2_dist_weights_lnempl1=lm(ln_empl_po ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru, weights=pol_02_au_ru$dist_weights)
summary(lm2_dist_weights_lnempl1)

lm2_dist_weights_lnempl2=lm(ln_empl_po ~ RU + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + gmina_wiej + gmina_mw + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru, weights=pol_02_au_ru$dist_weights)
summary(lm2_dist_weights_lnempl2)



### EMPL. / POP. (LOG.) - LAT. / LONG.

lm2_latlong_weights_lnempl1=lm(ln_empl_po ~ RU + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru, weights=pol_02_au_ru$dist_weights)
summary(lm2_latlong_weights_lnempl1)

lm2_latlong_weights_lnempl2=lm(ln_empl_po ~ RU + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + gmina_wiej + gmina_mw + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru, weights=pol_02_au_ru$dist_weights)
summary(lm2_latlong_weights_lnempl2)



### APP. / JOB (LOG.) - AIR DISTANCE

lm2_dist_weights_lnapp1=lm(ln_app_job ~ RU + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru, weights=pol_02_au_ru$dist_weights)
summary(lm2_dist_weights_lnapp1)

lm2_dist_weights_lnapp2=lm(ln_app_job ~ RU + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru, weights=pol_02_au_ru$dist_weights)
summary(lm2_dist_weights_lnapp2)



### APP. / JOB (LOG.) - LAT. / LONG.

lm2_latlong_weights_lnapp1=lm(ln_app_job ~ RU + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru, weights=pol_02_au_ru$dist_weights)
summary(lm2_latlong_weights_lnapp1)

lm2_latlong_weights_lnapp2=lm(ln_app_job ~ RU + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_au_ru + RU*dist_au_ru, data=pol_02_au_ru, weights=pol_02_au_ru$dist_weights)
summary(lm2_latlong_weights_lnapp2)



### ADVERT. - AIR DISTANCE

lm2_dist_weights_advt1_qp=glm(advert ~ RU + dist_au_ru + RU*dist_au_ru, family="quasipoisson", data=pol_02_au_ru, weights=pol_02_au_ru$dist_weights)
summary(lm2_dist_weights_advt1_qp)

lm2_dist_weights_advt2_qp=glm(advert ~ RU + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + dist_au_ru + RU*dist_au_ru, family="quasipoisson", data=pol_02_au_ru, weights=pol_02_au_ru$dist_weights)
summary(lm2_dist_weights_advt2_qp)



### ADVERT. - LAT. / LONG.

lm2_latlong_weights_advt1_qp=glm(advert ~ RU + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_au_ru + RU*dist_au_ru, family="quasipoisson", data=pol_02_au_ru, weights=pol_02_au_ru$dist_weights)
summary(lm2_latlong_weights_advt1_qp)

lm2_latlong_weights_advt2_qp=glm(advert ~ RU + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_au_ru + RU*dist_au_ru, family="quasipoisson", data=pol_02_au_ru, weights=pol_02_au_ru$dist_weights)
summary(lm2_latlong_weights_advt2_qp)



###########################
### PRUSSIA / AUSTRIA COMPARISON
###########################

### EMPL. / POP. (LOG.) - AIR DISTANCE

lm3_dist_weights_lnempl1=lm(ln_empl_po ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au, weights=pol_03_pr_au$dist_weights)
summary(lm3_dist_weights_lnempl1)

lm3_dist_weights_lnempl2=lm(ln_empl_po ~ AU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + gmina_wiej + gmina_mw + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au, weights=pol_03_pr_au$dist_weights)
summary(lm3_dist_weights_lnempl2)



### EMPL. / POP. (LOG.) - LAT. / LONG.

lm3_latlong_weights_lnempl1=lm(ln_empl_po ~ AU + GER1918 + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au, weights=pol_03_pr_au$dist_weights)
summary(lm3_latlong_weights_lnempl1)

lm3_latlong_weights_lnempl2=lm(ln_empl_po ~ AU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + gmina_wiej + gmina_mw + academ_app + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au, weights=pol_03_pr_au$dist_weights)
summary(lm3_latlong_weights_lnempl2)



### APP. / JOB (LOG.) - AIR DISTANCE

lm3_dist_weights_lnapp1=lm(ln_app_job ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au, weights=pol_03_pr_au$dist_weights)
summary(lm3_dist_weights_lnapp1)

lm3_dist_weights_lnapp2=lm(ln_app_job ~ AU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au, weights=pol_03_pr_au$dist_weights)
summary(lm3_dist_weights_lnapp2)



### APP. / JOB (LOG.) - LAT. / LONG.

lm3_latlong_weights_lnapp1=lm(ln_app_job ~ AU + GER1918 + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au, weights=pol_03_pr_au$dist_weights)
summary(lm3_latlong_weights_lnapp1)

lm3_latlong_weights_lnapp2=lm(ln_app_job ~ AU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_au + AU*dist_pr_au, data=pol_03_pr_au, weights=pol_03_pr_au$dist_weights)
summary(lm3_latlong_weights_lnapp2)



### ADVERT. - AIR DISTANCE

lm3_dist_weights_advt1_qp=glm(advert ~ AU + GER1918 + dist_pr_au + AU*dist_pr_au, family="quasipoisson", data=pol_03_pr_au, weights=pol_03_pr_au$dist_weights)
summary(lm3_dist_weights_advt1_qp)

lm3_dist_weights_advt2_qp=glm(advert ~ AU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + dist_pr_au + AU*dist_pr_au, family="quasipoisson", data=pol_03_pr_au, weights=pol_03_pr_au$dist_weights)
summary(lm3_dist_weights_advt2_qp)



### ADVERT. - LAT. / LONG.

lm3_latlong_weights_advt1_qp=glm(advert ~ AU + GER1918 + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_au + AU*dist_pr_au, family="quasipoisson", data=pol_03_pr_au, weights=pol_03_pr_au$dist_weights)
summary(lm3_latlong_weights_advt1_qp)

lm3_latlong_weights_advt2_qp=glm(advert ~ AU + GER1918 + ln_revenue + ln_pop_den + city_powia + migr_avg + unempl_avg + academ_app + lnpop + LocationLa + LocationLo + LocLaSq + LocLoSq + LocationLa:LocationLo + LocLaSq:LocationLo + LocationLa:LocLoSq + LocLaCu + LocLoCu + dist_pr_au + AU*dist_pr_au, family="quasipoisson", data=pol_03_pr_au, weights=pol_03_pr_au$dist_weights)
summary(lm3_latlong_weights_advt2_qp)



###########################
### MATCHING: ADDITIONAL INFORMATION
###########################

###########################
### GRAPHS OF THE PRUSSIA / RUSSIA COMPARISON
###########################

par(mfrow= c(3, 1))
plot(matching1_dv1, type = "jitter")
plot(matching1_dv2, type = "jitter")
plot(matching1_dv3, type = "jitter")

plot(matching1_dv1, type = "hist")
plot(matching1_dv2, type = "hist")
plot(matching1_dv3, type = "hist")



###########################
### GRAPHS OF THE AUSTRIA / RUSSIA COMPARISON
###########################

par(mfrow= c(3, 1))
plot(matching2_dv1, type = "jitter")
plot(matching2_dv2, type = "jitter")
plot(matching2_dv3, type = "jitter")

plot(matching2_dv1, type = "hist")
plot(matching2_dv2, type = "hist")
plot(matching2_dv3, type = "hist")



###########################
### GRAPHS OF THE PRUSSIA / AUSTRIA COMPARISON
###########################

par(mfrow= c(3, 1))
plot(matching3_dv1, type = "jitter")
plot(matching3_dv2, type = "jitter")
plot(matching3_dv3, type = "jitter")

plot(matching3_dv1, type = "hist")
plot(matching3_dv2, type = "hist")
plot(matching3_dv3, type = "hist")


